package Ascore_dev;

require Exporter;

our $VERSION = '1.00';
our(@ISA, @EXPORT);
@ISA = qw(Exporter);
@EXPORT = qw(ascore convert_sty);

use strict;
use warnings;
use Math::Combinatorics;
use Math::NumberCruncher;
use Data::Dumper;


sub ascore{
	my ($peptide, $parent_mass, $z, $masses,
		$intensities, $precision, $limit_mz_sup, $limit_mz_inf)= @_;
	
	#print "peptide standard notation: $peptide\n 
	#print "parent: $parent_mass \n charge: $z \n $masses \n $intensities \n";
	
	#El numero de sitios de fosforilacion ha de ser mayor que el de 
	#fosforilaciones para que el Ascore tenga sentido.
	#Si el Ascore no tiene sentido, lo devuelvo como 1000
	my $num_sitios_fosforilacion=0;
	my $num_sitios_deshidratacion=0;
	my $num_sitios_modificacion=0;
	my $num_fosforilaciones=0;
	my $num_deshidrataciones=0;
	my $num_modificaciones=0;
	
	my $peptido_conteo_modificaciones = $peptide;
	while($peptido_conteo_modificaciones=~s/(\(21\)|\(21\,|21\))/x/){
		$num_fosforilaciones++;
	}
	while($peptido_conteo_modificaciones=~s/(\(23\)|\(23\,|23\))/x/){
		$num_deshidrataciones++;
	}
	while($peptido_conteo_modificaciones=~s/^Y/x/){
		$num_sitios_fosforilacion++;
		$num_sitios_modificacion++;
	}
	while($peptido_conteo_modificaciones=~s/^[ST]/x/){
		$num_sitios_fosforilacion++;
		$num_sitios_deshidratacion++;
		$num_sitios_modificacion++;
	}
	while($peptido_conteo_modificaciones=~s/Y(\w)/x$1/){
		$num_sitios_fosforilacion++;
		$num_sitios_modificacion++;
	}
	while($peptido_conteo_modificaciones=~s/[ST](\w)/x$1/){
		$num_sitios_fosforilacion++;
		$num_sitios_deshidratacion++;
		$num_sitios_modificacion++;
	}
	
	my %result;
	#Si el numero de sitios de modificacion = numero de modificaciones
	$num_modificaciones = $num_fosforilaciones + $num_deshidrataciones;
	if($num_sitios_modificacion==$num_modificaciones){
		$result{'phosphorylations'} = $num_fosforilaciones;
		$result{'dehydratations'} = $num_deshidrataciones;
		$result{'peptide'} = $peptide;
		if($num_modificaciones==1){
			$result{'ascore'} = [1000];
		}
		elsif($num_modificaciones==2){
			$result{'ascore'} = [1000,999];
		}
		elsif($num_modificaciones==3){
			$result{'ascore'} = [1000,999,998];
		}
		return %result;
	}
	#Caso especial: numero de deshidrataciones = numero posibles sitios de deshidratacion
	if($num_deshidrataciones>0){
		if($num_sitios_deshidratacion == $num_deshidrataciones){
			if($num_modificaciones==1){
				$result{'phosphorylations'} = $num_fosforilaciones;
				$result{'dehydratations'} = $num_deshidrataciones;
				$result{'peptide'} = $peptide;
				$result{'ascore'} = [1000];
				return %result;
			}
			else{
				#Caso no resuelto: num.DHs == num.sitios.DH Y num_mods>=2
				$result{'phosphorylations'} = $num_fosforilaciones;
				$result{'dehydratations'} = $num_deshidrataciones;
				$result{'peptide'}=$peptide;
				if($num_modificaciones==2){
					$result{'ascore'} = [0,0];
				}
				elsif($num_modificaciones==3){
					$result{'ascore'} = [0,0,0];
				}
				print "CASO NO RESUELTO \n";
			}
		}
	}

	if($num_modificaciones>0){
		%result = process_ascore($peptide, $precision, $masses,
								 $intensities, $z, $limit_mz_sup, $limit_mz_inf);
	}
	else{
		$result{'peptide'} = $peptide;	
	}
	$result{'phosphorylations'} = $num_fosforilaciones;
	$result{'dehydratations'} = $num_deshidrataciones;
	return %result;
}


sub process_ascore{
	my $peptido = $_[0];		#en formato [styuv]
	my $precision = $_[1];		#en unidades m/z
	my $masses = $_[2];			#referencia a array de masas
	my $intensities = $_[3];	#referencia a array intensidades
	my $z = $_[4];				#carga
	my $limit_mz_sup = $_[5];			
	my $limit_mz_inf = $_[6];
	my @masas = split /\|/,$masses;
	my @intensidades = split /\|/, $intensities;
	my @peptidos_posibles = peptide_collection($peptido); 
	
	## VARIABLES USADAS
	# @pepts_by_score <- secuencias en formato original ordenadas por score de + a -
	# $best_pept <-----  secuencia peptido de mayor ascore 
	# @complements  <--  $site;
	# @col_peptidos <--  [@peptide_aas];   AoA de "peptido trozeados" formato stuv
	# @col_psites   <--  [@temp];     	   AoA de sitios de fosforilacion (incluido pept referencia)
	# @psite_pairs  <--  parejas compltarias para cada psite guardadas como indices de col_psites;
	# @best_psites                         sitios de fosfo en pept de mayor score
	# @levels_max_diff	<---  levels (10) de maxima diferencia de scores entre el primero y todos
	# $num_psites = $#best_psites+1;
	# $num_secuencias = $#col_psites+1;	
		
	#Obtencion de los Peptide Scores
	my %peptide_scores;
	my %matches_ascore;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos = %{prever_fragmentos($pep, $z, $limit_mz_sup, $limit_mz_inf)};
		$matches_ascore{$pep} = matches_ascore(\%fragmentos_teoricos, \@masas, 
											   \@intensidades, $precision, $pep);
	}
	
	my %weigths = (1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25, 10=>0.25);
	my $suma_niveles = 7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	
	foreach my $pep(keys %matches_ascore){
		my $suma_datos_matches = 0;
		my $score = 0;
		for (my $lvel=1; $lvel<11; $lvel++){
			$score += $weigths{$lvel} * $matches_ascore{$pep}{$lvel}{score};
			$matches_ascore{$pep}{$lvel}{score_ponderado} = ($weigths{$lvel}*$matches_ascore{$pep}{$lvel}{score}/$suma_niveles);
		}
		$peptide_scores{$pep} = $score/$suma_niveles;
	}
	#print Dumper %matches_ascore;
	
	#Ordeno peptidos segun Peptide score (de mayor a menor)
	my @pepts_by_score = sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $best_pept = $pepts_by_score[0];
	
	#print (my $so = $pepts_by_score[0]);
	#print "\n", $peptide_scores{$so}, "\n" ;
	#print (my $s1 = $pepts_by_score[-1]);
	#print "\n", $peptide_scores{$s1}, "\n";
	
	#Convierto peptidos en formato [styuv]
	my @col_peptidos;
	my @peptide_aas;
	foreach my $pept (@pepts_by_score){
		@peptide_aas = split //, convert_sty($pept);
		push @col_peptidos, [@peptide_aas];
	}

	#Detecto puntos de modificacion	
	my @col_psites;
	foreach my $peptide (@col_peptidos){
		my @temp;
		for(my $i=0; $i<$#{$peptide}+1; $i++){
			if(@{$peptide}[$i]=~/[styuv]/){
				push @temp, $i;
            }
		}
		push @col_psites, [@temp];
	}
	
	#Encuentro Peptidos-2 para cada modificacion
	#Los peptidos-2 se escogen con este criterio: se busca un peptido con el mayor 
	# peptide-score posible y que mantiene todos los puntos de forforilacion
	# del peptido de referencia excepto la fosforilacion estudiada
	#Otra forma de verlo, se buscan peptidos con ese punto cambiado 
	# pero que el numero de coincidencias sea = total_psites-1
	
	my @best_psites = @{$col_psites[0]};
	my $num_psites = $#best_psites + 1;
	my $num_secuencias = $#col_psites + 1;
	my @psite_pairs;
	
	#print "-\n";
	#print "\nref es : ", Dumper @{$best_psites};
	#print "pnum = ", $#{$best_psites}+1, "\n";
		
	my $psites_of_i;									# puntos de fosforilacion de la secuencia i
	for (my $site=0; $site<$num_psites; $site++){   	#para cada modificacion en el peptido de referencia
		for (my $i=1; $i<$num_secuencias; $i++){		#para cada una de las posibles secuencias
			$psites_of_i = $col_psites[$i];				#obten sus psites
			#mira si el sitio de referencia no esta fosforilado en la posible pareja
			if (!($best_psites[$site] ~~ @{$psites_of_i})){		
				my $count = 0;							#si es asi,
				foreach my $my_psite (@best_psites){	#calcula cuantos coinciden...
					if($my_psite ~~ @{$psites_of_i}){		
						$count += 1;			
					}
				}
				if ($count == $num_psites-1){			#si coinciden todos menos 1
					push @psite_pairs, $i;    		    #guarda el idx de la secuencia pareja
					#print "for site $site ";
					#print "la pareja es ", $i , "\n";
					last;								#y pasa a otro punto de fosforilacion
				}
			} 
		}
	}
	
	#Obtengo niveles maxima diferencia entre Peptide Scores
	#(P1 respecto a peptido2_f1, peptido2_f2 y peptido2_f3)
	my @levels_max_diff;
	my %datos_bp = %{$matches_ascore{$best_pept}};
	
	foreach my $mypeptide (@pepts_by_score){
		my $diferencia = 0;
		my $level;
		for(my $i=1; $i<11; $i++){
			my $dif = $datos_bp{$i}{score} - ${$matches_ascore{$mypeptide}}{$i}{score};
			if($dif > $diferencia){
				$diferencia = $dif;
				$level = $i;
			}
		}
		if($diferencia==0){		#Cuando el pep1 y el pep2 son identicos;
			$level = 5;			# uso el 5 de forma arbitraria 
		}
		#print "mayor diferencia es $diferencia, a nivel $level\n";
		push @levels_max_diff, $level;
	}	
		
	#Encuentro puntos complementarios a las f1,f2 y f3 en cada peptido
	#Por punto complementario se entiende al aminoacido que este modificado en p2 pero NO en p1
	my @complements;
	for(my $i=0; $i < $num_psites; $i++){           		   #de psite en psite del peptido referencia
		my @couple_sites = @{$col_psites[$psite_pairs[$i]]};   #cojo los sites del peptido pareja	
		foreach my $site (@couple_sites){					   #busco un site en el peptido pareja correspondiente
			#print "site $site para index $i\n";
			if (!($site ~~ @best_psites) and ($site != $best_psites[$i])){		#que no este en la referencia
				#print "para site $i el complementario es $site\n";
				push @complements, $site;
				last;
			}
		}
	}
	
	#print "ln 268 col_psites  --->\n ", Dumper @col_psites;
	#print "\n";
	#print "numero sitios fosforilados ", $num_psites, "\n";
	#print "numero de secuencias $num_secuencias\n";
	#print "site pairs\n", Dumper @psite_pairs;
	#print "ref psites\n", Dumper @best_psites;
	#print "\n";
	#print "complementos: \n", Dumper @complements;
	#print "\n";
	#return;
	
	#Calculo fragmentos	
	my @fingerprints_list;
	my $peptide_compl;
	my %fingerp_ref =  %{prever_fragmentos($pepts_by_score[0], $z,
										   $limit_mz_sup, $limit_mz_inf)};
	#print "referencia ", Dumper %fingerp_ref;
	
	my @frag_ref;
	my @frag_compl;
	my %fingerp_compl;
	my $start;
	my $end;
	my $idx = 0;
	foreach my $compsite (@complements){
		#print "compsite $compsite\n";
		
		if ($best_psites[$idx] > $compsite){
			$start = $compsite + 1;
			$end = $best_psites[$idx] + 1;
		}
		else {
			$start = $best_psites[$idx] + 1;
			$end = $compsite + 1;
		}
		
		$peptide_compl = $pepts_by_score[$psite_pairs[$idx]];
		
		print "peptido referencia ---> $pepts_by_score[0]\n";
		print "peptido_compl $idx ---> $peptide_compl\n\n";
		print "start-end, $start, $end\n";
		
		%fingerp_compl =  %{prever_fragmentos($peptide_compl, $z, $limit_mz_sup, $limit_mz_inf)};

		push @frag_ref, \%{get_specific_frags(convert_sty($pepts_by_score[0]),
						 					    \%fingerp_ref, $start, $end)};

		push @frag_compl, \%{get_specific_frags(convert_sty($peptide_compl),
												\%fingerp_compl, $start, $end)};
		$idx += 1;
	}
		
	#Obtengo datos-ascore: con matches de los fragmentos site-determining
	my @ascores;
	for (my $psite=0; $psite<$num_psites; $psite++){   
		my $peptide_compl = $pepts_by_score[$psite_pairs[$psite]];
		my $level = $levels_max_diff[$psite_pairs[$psite]];
		#
		my %dat_ascore_bp = %{matches_ascore($frag_ref[$psite], \@masas,
											  \@intensidades, $precision, $best_pept)};
		my %dat_ascore_compl = %{matches_ascore($frag_compl[$psite], \@masas,
												\@intensidades, $precision, $peptide_compl)};
		#
		my $matches_ref = $#{$dat_ascore_bp{$level}{tags_match_intervalo}} + 1;
		my $matches_compl = $#{$dat_ascore_compl{$level}{tags_match_intervalo}} + 1;
		#
		my @array = keys(%{$frag_ref[$psite]});
		my $size = 1 + $#array;
		my $proporcion = $matches_ref."//".$size." (".$matches_compl.")";
		
		
		my $value = $dat_ascore_bp{$level}{score} -
				     $dat_ascore_compl{$level}{score};
		
		if ($value < 0.1){$value = 0;}
		push @ascores, $value;
		
		#print "level_max, $level\n";
	}
	#print  "mi ascore =", Dumper @ascores;
	#print "\n";
	
	my %result;
	$result{'peptide'} = $best_pept;
	$result{'peptide_score'} = $peptide_scores{$best_pept};
	$result{'ascore'} = [@ascores];
	
	return %result;	
}

sub peptide_collection{
	my $pep = $_[0];
	#Asumo dos cosas:
	#	1/Solo hay un chemtag por cada aminoacido
	#	2/Los chemtag no cambian de sitio
	my @posibles_puntos_fosforilacion;
	my $num_fosforilaciones = 0;
	my @posibles_puntos_DH;
	my $num_DH = 0;
	my $num_chemtag = 0;

	my @aas;
	my @mods;
	my $aa_num = 0;

	while($pep=~/\w/){
		my $aa;
		my $mod;
		if($pep=~s/^(\w)\((\w+)\)//){
			$aa = $1;
			$mod = $2;
			push @aas, $aa;
			push @mods, $mod;
			#Calculo el numero de fosfos/chemtags/DH
			if($mod==21){
				$num_fosforilaciones++;
			}
			elsif($mod==23){
				$num_DH++;
			}
			elsif($mod==214 or $mod==737){
				$num_chemtag++;
				if($aa_num>0){#Si hay chemtag en un aa no N-terminal, no se puede modificar con nada mas
					$aa_num++;
					next;
				}
			}
			#Recopilo los posibles puntos de fosfo/DH
			if($aa=~/[S,T,Y]/){
				push @posibles_puntos_fosforilacion,$aa_num;
			}
			if($aa=~/[S,T]/){
				push @posibles_puntos_DH,$aa_num;
			}
		}
		else{
			$pep=~s/^(\w)//;
			$aa=$1;
			$mod='';
			push @aas,$aa;
			push @mods,$mod;
			#Recopilo los posibles puntos de fosfo/chemtag/DH
			if($aa=~/[S,T,Y]/){
				push @posibles_puntos_fosforilacion,$aa_num;
			}
			if($aa=~/[S,T]/){
				push @posibles_puntos_DH,$aa_num;
			}
		}
		$aa_num++;
	}
		
	#Obtengo las diferentes combinaciones de fosforilados, deshidratados y chemtag
	my @combinaciones_fosfo;
	my @combinaciones_DH;
	my @combinaciones_chemtag;
	
	my $combinat1 = Math::Combinatorics->new(count => $num_fosforilaciones,data => [@posibles_puntos_fosforilacion],);
	while(my @combo1 = $combinat1->next_combination){
		push @combinaciones_fosfo,\@combo1;
	}
	my $combinat2 = Math::Combinatorics->new(count => $num_DH,data => [@posibles_puntos_DH],);
		while(my @combo2 = $combinat2->next_combination){
			push @combinaciones_DH,\@combo2;
	}

	my @combinaciones_peptidos;
	if(@combinaciones_fosfo and @combinaciones_DH){
		#print "Hay phospho, DH y chemtag\t";
		foreach my $fos(@combinaciones_fosfo){
			my %usados;
			my @fos=@$fos;
			foreach my $f(@fos){
				$usados{$f}=1;
			}
			LABEL1:
			foreach my $dh(@combinaciones_DH){
				my @dhs=@$dh;
				foreach my $d(@dhs){
					if(defined $usados{$d}){
						next LABEL1;
					}
				}
				#Genero un nuevo array de modificaciones para construir el peptido modificado
				my @mods_new;
				for(my $i=0;$i<$#mods+1;$i++){#Recupero las modificaciones no phospho o DH (chemtag,MetOx,...)
					if($mods[$i]){
						unless(($mods[$i]==21)or($mods[$i]==23)){
							$mods_new[$i]='('.$mods[$i].')';
						}
					}
				}
				foreach(@fos){
					if($mods_new[$_]){
						my $t=$mods_new[$_];
						$t=~s/[\(,\)]//g;
						$mods_new[$_]='('.$t.',21)';	
					}
					else{
						$mods_new[$_]='(21)';
					}
				}
				foreach(@dhs){
					if($mods_new[$_]){
						my $t=$mods_new[$_];
						$t=~s/[\(,\)]//g;
						$mods_new[$_]='('.$t.',23)';	
					}
					else{
						$mods_new[$_]='(23)';
					}
				}
				#Genero el peptido
				my $peptide_string = '';
				for(my $i=0; $i<$#aas+1; $i++){
					if($mods_new[$i]){
						$peptide_string=$peptide_string.$aas[$i].$mods_new[$i];
					}
					else{
						$peptide_string=$peptide_string.$aas[$i];
					}
				}
				push @combinaciones_peptidos,$peptide_string;
			}
		}
	}
	elsif(@combinaciones_fosfo){
		#print "solo fosfo\t";
		foreach my $fos(@combinaciones_fosfo){
			my @fos=@$fos;
			# Nuevo array de modificaciones para construir el peptido modificado
			my @mods_new;
			# Recupero las modificaciones no phospho o DH (chemtag,MetOx,...)
			for(my $i=0;$i<$#mods+1;$i++){
				if($mods[$i]){
					unless(($mods[$i]==21)or($mods[$i]==23)){
						$mods_new[$i] = '('.$mods[$i].')';
					}
				}
			}
			foreach(@fos){
				if($mods_new[$_]){
					my $t=$mods_new[$_];
					$t=~s/[\(,\)]//g;
					$mods_new[$_] = '('.$t.',21)';	
				}
				else{
					$mods_new[$_] = '(21)';
				}
			}
			#Genero el peptido
			my $peptide_string = '';
			for(my $i=0; $i<$#aas+1; $i++){
				if($mods_new[$i]){
					$peptide_string = $peptide_string.$aas[$i].$mods_new[$i];
				}
				else{
					$peptide_string = $peptide_string.$aas[$i];
				}
			}
			push @combinaciones_peptidos,$peptide_string;
		}
	}
	elsif(@combinaciones_DH){
		#print "solo DH\t";
		foreach my $dhs(@combinaciones_DH){
			my @dhs=@$dhs;
			#Genero nuevo array de modificaciones para construir el peptido modificado
			my @mods_new;
			#Recupero las modificaciones no phospho o DH (chemtag,MetOx,...)
			for(my $i=0; $i<$#mods+1; $i++){
				if($mods[$i]){
					unless(($mods[$i]==21) or ($mods[$i]==23)){
						$mods_new[$i] = '('.$mods[$i].')';
					}
				}
			}
			foreach(@dhs){
				if($mods_new[$_]){
					my $t=$mods_new[$_];
					$t=~s/[\(,\)]//g;
					$mods_new[$_]='('.$t.',23)';
				}
				else{
					$mods_new[$_]='(23)';
				}
			}
			#Genero el peptido
			my $peptide_string='';
			for(my $i=0;$i<$#aas+1;$i++){
				if($mods_new[$i]){
					$peptide_string = $peptide_string.$aas[$i].$mods_new[$i];
				}
				else{
					$peptide_string = $peptide_string.$aas[$i];
				}
			}
			push @combinaciones_peptidos, $peptide_string;
		}
	}
	else{
		push @combinaciones_peptidos,$_[0];
	}
	return @combinaciones_peptidos;	
}

sub prever_fragmentos{
	# Predict a spectra using a sequence. Uses the 20 aa and allows phosphorilations
	#on STY (+80:s,t,y), water loss on ST (-18:u,v), and methionine oxidation (+18:m).
	#Produces a hash containing the {label}->mass.
	#It calculates the y and b series for the charge state, 
	#ARGUMENTS: there are 3 args:
	#       - $_[0]-> sequence
	#       - $_[1]-> charge state (1,2 or 3)
	#       - $_[2]-> mz upper limit for fragment generation
	#RETURN: a hash with {label}->mass.
	my $peptide = $_[0];
	my $z = $_[1];
	my $limit_mz_sup = $_[2];
	my $limit_mz_inf = $_[3];
	my %masas_aa = (
	A=>71.03711,  V=>99.06841,  L=>113.08406, I=>113.08406, P=>97.05276,
	M=>131.0405,  F=>147.0684,  W=>186.07932, G=>57.02146,  S=>87.03203,
	T=>101.04768, C=>160.03069, Y=>163.06332, N=>114.04293, Q=>128.05858,
	D=>115.02694, E=>129.04259, H=>137.05891, K=>128.09496, R=>156.1011);
	#la C esta CAM!, la considero PTM fija
	
	my @aas;
	my @mods;
	my @masas_aas;

	while($peptide=~/\w/){
		my $aa;
		my $mod;
		if($peptide=~s/^(\w)\((\w+\,\w+)\)//){
			$aa = $1;
			$mod = $2;
			push @aas, $aa;
			push @mods, $mod;
		}
		elsif($peptide=~s/^(\w)\((\w+)\)//){
			$aa = $1;
			$mod = $2;
			push @aas, $aa;
			push @mods, $mod;
		}
		else{
			$peptide=~s/^(\w)//;
			$aa = $1;
			$mod = '';
			push @aas, $aa;
			push @mods, $mod;
		}
	}
	my %masas_mod = (21 => 79.9799,         #fosfato
                     23 => -18.010565,      #dh
                     214=> 144.102063,      #itraq
                     737=> 229.162932,      #tmt
                     35 => 15.994915);      #oxidacion
     
    for(my $i=0; $i<$#aas+1; $i++){
		if($mods[$i]){
			if($mods[$i]=~/\,/){
				my ($m1, $m2) = split /\,/, $mods[$i];
				$masas_aas[$i] = $masas_mod{$m1} + $masas_mod{$m2} + $masas_aa{$aas[$i]};
			}
			else{
				$masas_aas[$i] = $masas_mod{$mods[$i]} + $masas_aa{$aas[$i]};
				#print $mods[$i],"\n";
			}
		}
		else{
			$masas_aas[$i] = $masas_aa{$aas[$i]};
		}
	}
	
	#La carga maxima permitida es 1 para z=1 o 2, 
	my $fragment_max_charge;
	if($z<3){
		$fragment_max_charge = 1;
	}
	else{
		$fragment_max_charge = $z-1;
	}

	my $acumulado = 0;
	my (@serie_y_1, @serie_b_1, @serie_y_2, @serie_b_2, @serie_y_3, @serie_b_3);

	for (my $i=0; $i<$#aas+1; $i++){
		$serie_b_1[$i] = $masas_aas[$i] + $acumulado+1;
		$acumulado += $masas_aas[$i];
		}
	my $M = $acumulado + 19;						####ojo este 19!!!!!!
	for(my $i=0;$i<$#aas;$i++){
		$serie_y_1[$i] = $M - $serie_b_1[($#aas-$i)] + 1;
		}
	$serie_y_1[$#aas] = $M - $serie_b_1[0] + 1;
	if($_[1]>1){
		for (my $i=0; $i<$#aas+1; $i++){
			$serie_y_2[$i] = ($serie_y_1[$i]+1)/2;
			if (defined $serie_b_1[$i]){
				$serie_b_2[$i] = ($serie_b_1[$i]+1)/2;
				}
			}
		if($_[1]>2){
			for (my $i=0; $i<$#aas+1; $i++){
				$serie_y_3[$i] = ($serie_y_1[$i]+1)/3;
				$serie_b_3[$i] = ($serie_b_1[$i]+1)/3;
			}		
		}	    	
	}

	my %fragmentos;
	for (my $i=0; $i<$#aas+1; $i++){
		my $j = $i+1; 
		if ($i<($#aas)){														# El b(n) bo existe
			my $num = $j;
			if($num<10){$num = '0'.$num};
			$fragmentos{"b$num+1"} = sprintf("%.3f", $serie_b_1[$i]);
			if(($_[1]>1)and($fragment_max_charge>1)){					# Si carga>1
				$fragmentos{"b$num+2"} = sprintf("%.3f", $serie_b_2[$i]);
				if(($_[1]>2) and ($fragment_max_charge>2)){				# Si carga>2
					$fragmentos{"b$num+3"} = sprintf("%.3f", $serie_b_3[$i]);
				}
			}	
		}
		if(($i>0)){	
			my $num = $i;
			if($num<10){$num='0'.$num};											# El y(0) bo existe
			$fragmentos{"y$num+1"} = sprintf("%.3f", $serie_y_1[$i]);
			 if(($_[1]>1)and($fragment_max_charge>1)){				# Si carga>1
				$fragmentos{"y$num+2"} = sprintf("%.3f", $serie_y_2[$i]);
				if(($_[1]>2)and($fragment_max_charge>2)){				# Si carga>2
					$fragmentos{"y$num+3"} = sprintf("%.3f", $serie_y_3[$i]);	
				}	
			}
		}	
	}
	#Solo se consideran los fragmentos entre limit_mz_sup e _inf
	foreach my $frag(keys %fragmentos){
		if(($fragmentos{$frag}>$limit_mz_sup)or($fragmentos{$frag}<$limit_mz_inf)){
			delete $fragmentos{$frag};
		}
	}
	return \%fragmentos;
}


sub binomial{
	#Binomial distribution
	#($n,$k,$p)->$k successes in $n tries, given a probability of $p	
	my $n = $_[0];
	my $k = $_[1];
	my $p = $_[2];
	my $binomial_acumulada=0;
	for(my $i=$k; $i<$n+1; $i++){
		my $binomial = Math::NumberCruncher::Binomial($n,$i,$p);
		$binomial_acumulada = $binomial_acumulada+$binomial;
	}
	return $binomial_acumulada;
}


sub convert_sty{
	my $pep=$_[0];
	#cambio fosfo o DH simples
	$pep=~s/S\(21\)/s/g;
	$pep=~s/T\(21\)/t/g;
	$pep=~s/Y\(21\)/y/g;
	$pep=~s/S\(23\)/u/g;
	$pep=~s/T\(23\)/v/g;
	#cambio fosfo o DH combinados con chemtag en Nt
	$pep=~s/S\((214|737)\,21\)/s/;
	$pep=~s/T\((214|737)\,21\)/t/;
	$pep=~s/Y\((214|737)\,21\)/y/;
	$pep=~s/S\(21\,(214|737)\)/s/;
	$pep=~s/T\(21\,(214|737)\)/t/;
	$pep=~s/Y\(21\,(214|737)\)/y/;
	$pep=~s/S\((214|737)\,23\)/u/;
	$pep=~s/T\((214|737)\,23\)/v/;
	$pep=~s/S\(23\,(214|737)\)/u/;
	$pep=~s/T\(23\,(214|737)\)/v/;
	#Elimino todas las demas modificaciones
	$pep=~s/\(\d+\)//g;
	return $pep;
}

sub get_specific_frags{
	#Genero un nuevo hash para poner los fragmentos que estan entre el primero y ultimo aa que se pueda fosforilar
	my $peptido = $_[0];
	my %fragmentos = %{$_[1]};
	my $primer_aa_fosforilable = $_[2];
	my $ultimo_aa_fosforilable = $_[3];
	my %fragmentos_filtrados;
	my $largo_pep = length $peptido;
	
	#print "peptido $peptido\n";
	#print "fragmentos ", Dumper $_[1];
	#print "largo $largo_pep\n";
	
	foreach my $frag (keys %fragmentos){
		my $frag2=(split /\+/,$frag)[0];
		if($frag2=~m/y/){
			$frag2=~s/y//;		#para las y, valen los fragmentos entre [largo-ultimo+1 , largo-primero]
			if(($frag2 >= ($largo_pep-$ultimo_aa_fosforilable+1))and
			   ($frag2 <= ($largo_pep-$primer_aa_fosforilable))){
				$fragmentos_filtrados{$frag} = $fragmentos{$frag};
			}
		}
		elsif($frag2=~m/b/){	#para las b, valen los fragmentos entre [primero , ultimo-1]
			$frag2=~s/b//;
			if(($frag2>=($primer_aa_fosforilable))and($frag2<=($ultimo_aa_fosforilable-1))){
				$fragmentos_filtrados{$frag}=$fragmentos{$frag};
			}
		}
	}
	%fragmentos = ();
	%fragmentos = %fragmentos_filtrados;	#Uso solo fragmentos en zona de interes		
	return \%fragmentos;
}


sub matches_ascore{
	my %fragmentos = %{$_[0]};
	my @masas= @{$_[1]};
	my @intensidades = @{$_[2]};
	my $fragment_error = $_[3];
	my $peptido = $_[4];
	
	#print "ascore_frag_error, $_[3]\n";
	#print "peptido, $peptido\n";
	
	my $num_fragmentos = keys %fragmentos;
	
	my $picos_x_sector;
	my %datos_num_picos;
	my @rangos_masas;
	my ($rango_masas_inf, $rango_masas_sup) = (200, 2000);
	for(my $i=$rango_masas_inf; $i<$rango_masas_sup; $i=$i+100){
		push @rangos_masas, $i;
	}
	for($picos_x_sector=1; $picos_x_sector<11; $picos_x_sector++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;
		my $correspondencias = 0;
		my $num_peaks_intervalo = 0;
		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior = $intervalo-100;
			for(my $i=0; $i<$#masas+1; $i++){
				if(($masas[$i]>$intervalo_inferior) and ($masas[$i]<=$intervalo)){
					$picos_intervalo{$masas[$i]} = $intensidades[$i];
					$num_peaks_intervalo++;
				}
			}
			
			my @masas_by_intensity;
			
			#ningun pico experimental en el intervalo
			if($num_peaks_intervalo == 0){
				next;
			}
			#un solo pico experimental en el intervalo
			elsif($num_peaks_intervalo==1){ 
				foreach my $key(keys %picos_intervalo ){
					push @masas_by_intensity, $key;
				}
			}
			#mas de 1 pico experimental en el intervalo
			else{	
				foreach my $mass (sort { $picos_intervalo{$b} <=> $picos_intervalo {$a}} keys %picos_intervalo ){
					push @masas_by_intensity, $mass;
				}
			}
			
			#Matching entre masas teoricas y experimentales
			my $masas_testadas = 0;
			foreach my $masas_exp(@masas_by_intensity){
				$masas_testadas++;
				foreach(keys %fragmentos){
					if(abs($fragmentos{$_}-$masas_exp) <= $fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo, $masas_exp;
						push @masas_teor_match_intervalo, $fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas == $picos_x_sector){
					last;
				}
			}
		}
		#$num_fragmentos=> trials
		#$correspondencias=>  successes
		#$prob=> numero picos por sector / 100

		my $probablity = ($picos_x_sector/100);
		my $binomial_acum_rango = binomial($num_fragmentos, $correspondencias, $probablity);
		my $score;
		if($binomial_acum_rango>0){
			$score = -10*(log($binomial_acum_rango)/log(10));
		}
		else{
			$score = 0;
		}
		#print "score, $score\n";
		$datos_num_picos{$picos_x_sector}{masas_exp_match_intervalo} = \@masas_exp_match_intervalo;
		$datos_num_picos{$picos_x_sector}{masas_teor_match_intervalo} = \@masas_teor_match_intervalo;
		$datos_num_picos{$picos_x_sector}{tags_match_intervalo} = \@tags_match_intervalo;
		$datos_num_picos{$picos_x_sector}{score} = $score;
	}
	return \%datos_num_picos;
}


return 1;

__END__

=head1 NAME

=head1 SYNOPSIS

=head1 DESCRIPTION

=head1 LIMITATIONS

=head1 SEE ALSO

=cut



